This computer lab addresses the implementation and analysis of reconstruction algorithms for compressive sensing. In particular, an application to Medical resonance imaging (MRI) is addressed. The goal is to show the connection between compressive sensing and denoising. To get an intuition, this Lab first explores synthetic 1D sparse signals
## sparse vector (nnz/length = 5/128) of class "dsparseVector"
## [1] . . . . . . . . . . . . . . . . . .
## [19] . 0.8 . . . . . . . . . . . . 0.2 . . .
## [37] . . . 0.6 . . . . . . . . 0.4 . . . . .
## [55] . . . . . . . . . . . . . . . . . .
## [73] . . . . . . . . . . . . . . . . . .
## [91] . . . . 1.0 . . . . . . . . . . . . .
## [109] . . . . . . . . . . . . . . . . . .
## [127] . .
data <- cbind.data.frame(y = as.numeric(x), x = 1:128)
ggplot(data, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2) +
xlab("position")+
ylab("x")
noise <- rnorm(128, mean = 0, sd = 0.05)
y <- x + noise
noisydata <- cbind.data.frame(y = as.numeric(y), x = 1:128)
ggplot(noisydata, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
ggtitle("Noisy y") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("position")+
ylab("x")
\[arg\min_{z \in \mathbb{R}^{n}} {1\over2}\left \| z - y \right \|_{2}^{2} + {1\over2}\lambda\left \| z \right \|_{2}^{2} \] Lets denote : \[ Lets \ denote : \\ \phi_{y}(u) = {1\over2}u^{2}\lambda + {1\over2}(u-y)^{2} \\ \phi_{y}^{'}(u^{*}) = {1\over2}u^{*}\lambda + {1\over2}(2u^{*}-2y) \\ \ = u^{*}(1 + \lambda) - y \\ \Rightarrow u^{*} = {1\over 1+\lambda}y \]
lambda <- c(0.01, 0.05, 0.1, 0.2)
xhat <- matrix(nrow = length(y), ncol = length(lambda))
y <- as.numeric(y)
for (i in 1:length(lambda)) {
for(j in 1:length(y)) {
xhat[j, i] <- (1/(1+lambda[i]))*y[j]
}
}
Lets plot xhat for = 0.1
denoisydata <- cbind.data.frame(y = xhat[,3], x = 1:128)
ggplot(denoisydata, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="green", fill=alpha("lightgreen", 0.3), alpha=0.7, shape=21, stroke=2) +
ggtitle("’l_2 norm denoising’") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("position")+
ylab("x")
Here we see that the solution is not sparse, this is due to the fact that wee use the l2 norm. Here the coefficient are shriken by {11 + } but they won’t go to zero.
\[ arg \min_{z \in \mathbb{R}^{n}} {1\over2}\left \| z - y \right \|_{2}^{2} + \lambda\left \| z \right \|_{1} = arg \min_{z \in \mathbb{R}^{n}}\phi_{y}(u) \]
\[\phi_{y}(u)\] is “separable” :
\[ \phi_{y}(u) = {1\over2} \sum_{i=1}^{n} (u_{i}-y_{i})^{2} +\lambda \sum_{i=1}^{n}(u_{i}) = \sum_{i=1}^{n} {1\over 2} (u_{i}-y_{i})^{2} + \lambda (u_{i}) = \sum_{i=1}^{n} \phi_{y_{i}}(u_{i}) \]
\[arg \min_{u_{1},u_{2},...,u_{n}}\phi_{y}(u) = \begin{pmatrix} arg \min_{u_{1}}\phi_{y_{1}}(u_{1}) \\ ... \\ arg \min_{u_{n}}\phi_{y_{n}}(u_{n}) \\ \end{pmatrix} = u^{*} \]
We need : \[arg \min_{u \in \mathbb{R}} {1\over2}(u - y )^{2} + \lambda\left \| u \right \| \]
first case \[u > 0\] : \[\phi_{y}(u) = \lambda u + {1\over2}(u-y)^{2}\] Lets try to find the minimum : \[\phi_{y}^{'}(u^{*}) = \lambda + (u-y) = 0 \\ \Leftrightarrow u_{1}^{*} = [y - \lambda]^{+} \] If \(u < 0\):
\[ \phi_{y}(u) = - \lambda u + \frac{1}{2t} + (u-y)^{2} \\ \phi_{y}^{'}(u) = 0 \\ u_{2}^{*} = [y+ \lambda t]^{-}\]
Equivalently, If $y $ :
\[y \in [- \lambda, \lambda ] \\ u_{1}^{*} = x + \lambda \\ u_{2}^{*} = 0\]
\[ \phi_{y}(u_{1}^{*}) = \lambda (y - \lambda) + \frac{1}{2} (y - \lambda - y)^{2} \\ \phi_{y}(u_{2}^{*}) = \frac{1}{2 \lambda} y^{2} \\ u_{1}^{*} = x - \lambda t \\ u_{2}^{*} = 0 \\ \phi_{y}(u_{1}^{*}) \leq \phi_{y}(u_{2}^{*}) \\ u^{*} = u_{1}^{*} \]
\[\hat{x} = SoftThresh(y, λ) = \begin{equation} \left\{ \begin{aligned} y_{l} + \lambda,& \ \ \ \text{if } y_{l} \leq \lambda\\ 0,& \ \ \ \text{if } \mid y_{l}\mid < \lambda\\ y_{l} - \lambda,& \ \ \ \text{if } y_{l} \geq \lambda\\ \end{aligned} \right. \end{equation}\]
The differece between Tichonov regularization with l1 norm and the Basis Pursuit (BP) algorithm is that Basis Pursuit has a fat A matrix behind the z in the first terrm of the function we wan’t to minimize. Basis Pursuit provides the sparse solution that Tikhonov regularization doesn’t.
SoftThresh <- function(y, lambda) {
xhat <- vector(length = length(y))
for( i in which((abs(y)<lambda), arr.ind = TRUE)) {
xhat[i] <- 0
}
for( i in which((y<= -lambda), arr.ind = TRUE) ) {
xhat[i] <- y[i] + lambda
}
for( i in which((y>=lambda), arr.ind = TRUE) ) {
xhat[i] <- y[i] - lambda
}
return(xhat)
}
(SoftThresh_y <- SoftThresh(y = -10:10, lambda = 2))
## [1] -8 -7 -6 -5 -4 -3 -2 -1 0 0 0 0 0 1 2 3 4 5 6 7 8
Lets Plot SoftThresh(u, λ) for u ∈ [−10, 10] and λ = 2
SoftThresh_data <- cbind.data.frame(y =SoftThresh_y , x = -10:10)
ggplot(SoftThresh_data, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="yellow", fill=alpha("lightyellow", 0.3), alpha=0.7, shape=21, stroke=2) +
ggtitle("’SoftThresh function, λ =num2str(lambda)]’") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("position")+
ylab("x")
When y is small compare to lambda, all the coefficients of will be 0. If y is big compare to lambda,the coefficients of will be pretty close to the coefficients of y.
SoftThresh_xhat <- matrix(nrow = length(y), ncol = length(lambda))
for (i in 1:length(lambda)) {
SoftThresh_xhat[, i] <- SoftThresh(y = y, lambda = lambda[i])
}
Lets plot xhat for = 0.1
SoftThresh_denoisy_data <- cbind.data.frame(y = SoftThresh_xhat[,3], x = 1:128)
ggplot(SoftThresh_denoisy_data, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="green", fill=alpha("lightgreen", 0.3), alpha=0.7, shape=21, stroke=2) +
ggtitle("’SoftThresh denoising’") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("position")+
ylab("x")
Here we can oberve a pretty good estimate but not perfect of the original vector x. If we take a look at the estimation we can see that indeed it is sparse but the hight of the non-zero coefficients is not the same as in the original vector.
We’ll now explore the connection between compressive sensing and denoising and the importance of choosing a good measurement matrix. To do so, we’ll observe the effect of regular and then random sampling of the signal in the frequency domain.
lets clear all variables :
rm(x, y, noise, data, noisydata, denoisydata, SoftThresh_data, SoftThresh_denoisy_data, xhat, SoftThresh_xhat)
(x <- sparseVector(x = ((1:5)/5), i = sample(128, 5), length=128))
## sparse vector (nnz/length = 5/128) of class "dsparseVector"
## [1] . . . . 1.0 . . . . . 0.2 0.8 . . . . . .
## [19] . . . . . . . . . . . . . . . . . .
## [37] . . . . . . . . . . . 0.6 . . . . . .
## [55] . . . . . . . . . . . . . . . . . .
## [73] . . . 0.4 . . . . . . . . . . . . . .
## [91] . . . . . . . . . . . . . . . . . .
## [109] . . . . . . . . . . . . . . . . . .
## [127] . .
data <- cbind.data.frame(y = as.numeric(x), x = 1:128)
ggplot(data, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2)+
xlab("position")+
ylab("x")
X <- fft(as.numeric(x))
n <- 128
Xu <- rep(0, n)
Xu[seq(from = 1, to = n, by = 4)] <- X[seq(from = 1, to = n, by = 4)]
xu <- fft(Xu, inverse=TRUE)*4
Let us denote T the matrix which corresponds to the discrete Fourier transform, i.e. X = T x. Recall that T is an orthogonal matrix. Let us denote TK the matrix which computes K out of n Fourier coefficients, i.e. XK = TKx, where XK is a K-length vector. Show that : \[x_{K} = T_{T}^{K}X_{K} = T_{T}^{K}T_{T}x_{K} \] is the orthogonal projection of x on the subspace \[I_{m}(T_{T}^{K})\]. Show that this is the minimum \[l_{2}\] norm solution. We knoz that T is an orthogonal matrix, so : \[TT^{t} = I\].
We know that :\[ X=Tx \Rightarrow T^{-1}X=T^{-1}Tx \Rightarrow x=T^{-1}X=x\]
realpart <- cbind.data.frame(y = Re(xu), x = 1:length(Re(xu)))
ggplot(realpart, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the real part", y = "Real part") +
theme(plot.title = element_text(hjust = 0.5))
imaginary_part <- cbind.data.frame(y = Im(xu), x = 1:length(Im(xu)))
ggplot(imaginary_part, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the imaginary part", y = "Imaginary part") +
theme(plot.title = element_text(hjust = 0.5))
So here we observe that xu semms to be periodic. We can see that the imaginary part is very close to zero so we can say that xu only has real part. The real part is periodic, it represent Xu periodicaly in the time space so we know the frequence. But we don’t know the shift of xu in the time space. This lack of information does not allow us to reconstruct the original signal from the result. The fact that we use equidistant points to generate Xu from X has periodise the signal xu.
Xr <- rep(0, n)
prm <- sample(128, 32)
Xr[prm] <- X[prm]
xr <- fft(Xr, inverse=TRUE) * 4
require(gridExtra)
## Loading required package: gridExtra
realpartr <- cbind.data.frame(y = Re(xr), x = 1:length(Re(xr)))
real_plotr <- ggplot(realpartr, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the real part", y = "Real part", x="position") +
theme(plot.title = element_text(hjust = 0.5))
imaginary_partr <- cbind.data.frame(y = Im(xr), x = 1:length(Im(xr)))
imaginary_plotr <- ggplot(imaginary_partr, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the imaginary part", y = "Imaginary part", x="position") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(real_plotr, imaginary_plotr, ncol=2)
by using random sampling we observe that the signal is not periodic anymore but still not very sparse we have now a signal denoising problem. The noise we see is actually incoherent aliasing that is contributed by the signal itself.
lets clear all variables :
rm(x, data, imaginary_part, realpart, imaginary_partr, realpartr, imaginary_plotr, real_plotr)
(x <- sparseVector(x = ((1:5)/5), i = sample(128, 5), length=128))
## sparse vector (nnz/length = 5/128) of class "dsparseVector"
## [1] . . . . . . . . . . . . . . . . . .
## [19] . . . . . . . . . . . 0.8 . . . . . .
## [37] . . . . 0.4 . . . . . . . . . . 0.6 . .
## [55] 1.0 . . . . . . . . . . . . . . . . .
## [73] . . . . . . . . . . . . . . . . 0.2 .
## [91] . . . . . . . . . . . . . . . . . .
## [109] . . . . . . . . . . . . . . . . . .
## [127] . .
data <- cbind.data.frame(y = as.numeric(x), x = 1:128)
ggplot(data, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2) +
xlab("position")+
ylab("x")
X <- fft(as.numeric(x))/sqrt(length(x))
Y <- rep(0, n)
prm <- sample(128, 32)
Y[prm] <- X[prm]
y <- fft(Y, inverse=TRUE)/sqrt(length(x))
require(gridExtra)
real_part <- cbind.data.frame(y = Re(y), t = 1:length(Re(y)))
real_plot <- ggplot(real_part, aes(x=t, y=y)) +
geom_segment( aes(x=t, xend=t, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the real part", y = "Real part", x="position") +
theme(plot.title = element_text(hjust = 0.5))
imaginary_part <- cbind.data.frame(y = Im(y), t = 1:length(Im(y)))
imaginary_plot <- ggplot(imaginary_part, aes(x=t, y=y)) +
geom_segment( aes(x=t, xend=t, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the imaginary part", y = "Imaginary part", x="position") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(real_plot, imaginary_plot, ncol=2)
SoftThreshComplex <- function(y, lambda) {
xhat <- vector(length = length(y))
for( i in which((Mod(y)<lambda), arr.ind = TRUE)) {
xhat[i] <- 0
}
for( i in which((Mod(y)>lambda), arr.ind = TRUE) ) {
xhat[i] <- y[i]*(Mod(y[i])-lambda) / Mod(y[i])
}
return(xhat)
}
Xhat <- Y
i <- 1
while(i < 101) {
Xhat <- Xhat
xhat <- fft(Xhat, inverse = TRUE)/sqrt(length(Xhat))
xhat <- SoftThreshComplex(y = xhat, lambda = 0.05)
Xhat <- (fft(xhat)/sqrt(length(xhat)))
Xhat <- Xhat*(Y==0) + Y
i <- i+1
}
real_part <- cbind.data.frame(y = Re(xhat), x = 1:length(Re(xhat)))
real_plot <- ggplot(real_part, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the real part", y = "Real part", x="position") +
theme(plot.title = element_text(hjust = 0.5))
imaginary_part <- cbind.data.frame(y = Im(xhat), x = 1:length(Im(xhat)))
imaginary_plot <- ggplot(imaginary_part, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the imaginary part", y = "Imaginary part", x="position") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(real_plot, imaginary_plot, ncol=2)
lambda <- c(0.01, 0.05, 0.1)
Xhat <- Y
xhat <- fft(Xhat, inverse = TRUE)/sqrt(length(xhat))
x <- as.numeric(x)
SoftThreshComplex_xhat <- matrix(nrow = length(xhat), ncol = length(lambda))
SoftThreshComplex_Xhat <- matrix(nrow = length(xhat), ncol = length(lambda))
error <- matrix(nrow = 100, ncol = length(lambda))
errorfun <- function(x, xhat){
return(sum(x^2-xhat^2))
}
for (i in 1:length(lambda)) {
SoftThreshComplex_Xhat[, i] <- Xhat
}
for (j in 1:100) {
SoftThreshComplex_Xhat <- SoftThreshComplex_Xhat
for (i in 1:length(lambda)) {
SoftThreshComplex_xhat[, i] <- fft(SoftThreshComplex_Xhat[,i], inverse = TRUE)/sqrt(length(xhat))
}
for (i in 1:length(lambda)) {
SoftThreshComplex_xhat[, i] <- SoftThreshComplex(y = SoftThreshComplex_xhat[, i], lambda = lambda[i])
}
for (i in 1:length(lambda)) {
SoftThreshComplex_Xhat[, i] <- fft(SoftThreshComplex_xhat[, i])/sqrt(length(xhat))
}
for (i in 1:length(lambda)) {
SoftThreshComplex_Xhat[, i] <- SoftThreshComplex_Xhat[, i]*(Y==0) + Y
error[j,i] <- errorfun(as.numeric(x), Re(SoftThreshComplex_xhat[, i]))
}
}
for (i in 1:length(lambda)) {
SoftThreshComplex_xhat[, i] <- fft(SoftThreshComplex_Xhat[,i], inverse = TRUE)/sqrt(length(xhat))
}
real_plot <- ggplot(data, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.3), alpha=0.7, shape=21, stroke=2)
real_part1 <- cbind.data.frame(y = Re(SoftThreshComplex_xhat[, 1]), x = 1:length(Re(SoftThreshComplex_xhat[, 1])))
real_plot1 <- ggplot(real_part1, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Estimation of x with lambda = 0.01", y = "Real part x_hat", x = "position") +
theme(plot.title = element_text(hjust = 0.5))
real_part2 <- cbind.data.frame(y = Re(SoftThreshComplex_xhat[, 2]), x = 1:length(Re(SoftThreshComplex_xhat[, 2])))
real_plot2 <- ggplot(real_part2, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Estimation of x with lambda = 0.05", y = "Real part x_hat", x = "position") +
theme(plot.title = element_text(hjust = 0.5))
real_part3 <- cbind.data.frame(y = Re(SoftThreshComplex_xhat[, 3]), x = 1:length(Re(SoftThreshComplex_xhat[, 3])))
real_plot3 <- ggplot(real_part3, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Estimation of x with lambda = 0.1", y = "Real part x_hat", x = "position") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(real_plot, real_plot1, real_plot2, real_plot3, ncol=1, nrow = 4)
We can observe above that the best aproximation of the original signal is obtain with lambda = 0.05. As lambda gets larger more coefficients are put to 0 and the approximation becomes more inaccurate because it can’t capture anymore all the important coefficients in order to reconstruct the original vector x. Conversely, as lambda gets smaller some noise are still present.
plot(error[,1], main = "Error for Lambda = 0.01", ylab = "Error", xlab = "Iterations")
plot(error[,2], main = "Error for Lambda = 0.05", ylab = "Error", xlab = "Iterations")
plot(error[,3], main = "Error for Lambda = 0.1", ylab = "Error", xlab = "Iterations")
From the mean square error plots we can observe that the lowest error is fro lambda = 0.05. As lambda gets larger the faster the error converges to 0 because a larger shrinkage is applied on the coefficients. We have to choose carefully the value of lambda and the number of iterations of the algorithm.
X <- fft(x)/sqrt(length(x))
Xu <- rep(0, 128)
Xu[seq(from = 1, to = 128, by = 4)] <- X[seq(from = 1, to = 128, by = 4)]
xu <- fft(Xu, inverse=TRUE)/sqrt(length(Xu))
realpart <- cbind.data.frame(y = Re(xu), x = 1:length(Re(xu)))
ggplot(realpart, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.3), alpha=0.7, shape=21, stroke=2) +
labs(title = "Plot of the real part", y = "Real part") +
theme(plot.title = element_text(hjust = 0.5))
Xhat <- Xu
xhat <- fft(Xhat, inverse = TRUE)/sqrt(length(Xhat))
x <- as.numeric(x)
SoftThreshComplex_xhat <- matrix(nrow = length(xhat), ncol = length(lambda))
SoftThreshComplex_Xhat <- matrix(nrow = length(xhat), ncol = length(lambda))
error <- matrix(nrow = 100, ncol = length(lambda))
errorfun <- function(x, xhat){
return(sum(x^2-xhat^2))
}
for (i in 1:length(lambda)) {
SoftThreshComplex_Xhat[, i] <- Xhat
}
for (j in 1:100) {
SoftThreshComplex_Xhat <- SoftThreshComplex_Xhat
for (i in 1:length(lambda)) {
SoftThreshComplex_xhat[, i] <- fft(SoftThreshComplex_Xhat[,i], inverse = TRUE)/sqrt(length(SoftThreshComplex_Xhat[,i]))
}
for (i in 1:length(lambda)) {
SoftThreshComplex_xhat[, i] <- SoftThreshComplex(y = SoftThreshComplex_xhat[, i], lambda = lambda[i])
}
for (i in 1:length(lambda)) {
SoftThreshComplex_Xhat[, i] <- fft(SoftThreshComplex_xhat[, i])/sqrt(length(SoftThreshComplex_xhat[, i]))
}
for (i in 1:length(lambda)) {
SoftThreshComplex_Xhat[, i] <- SoftThreshComplex_Xhat[, i]*(Xu==0) + Xu
error[j,i] <- errorfun(as.numeric(x), Re(SoftThreshComplex_xhat[, i]))
}
}
for (i in 1:length(lambda)) {
SoftThreshComplex_xhat[, i] <- fft(SoftThreshComplex_Xhat[,i], inverse = TRUE)/sqrt(length(xhat))
}
real_plot <- ggplot(data, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="red", fill=alpha("orange", 0.1), alpha=0.7, shape=21, stroke=2)+
labs(title = "true signal", y = "Real part of x", x = "position") +
theme(plot.title = element_text(hjust = 0.5))
real_part1 <- cbind.data.frame(y = Re(SoftThreshComplex_xhat[, 1]), x = 1:length(Re(SoftThreshComplex_xhat[, 1])))
real_plot1 <- ggplot(real_part1, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.1), alpha=0.7, shape=21, stroke=2) +
labs(title = "Estimation of x with lambda = 0.01", y = "Real part x_hat", x = "position") +
theme(plot.title = element_text(hjust = 0.5))
real_part2 <- cbind.data.frame(y = Re(SoftThreshComplex_xhat[, 2]), x = 1:length(Re(SoftThreshComplex_xhat[, 2])))
real_plot2 <- ggplot(real_part2, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.1), alpha=0.7, shape=21, stroke=2) +
labs(title = "Estimation of x with lambda = 0.05", y = "Real part x_hat", x = "position") +
theme(plot.title = element_text(hjust = 0.5))
real_part3 <- cbind.data.frame(y = Re(SoftThreshComplex_xhat[, 3]), x = 1:length(Re(SoftThreshComplex_xhat[, 3])))
real_plot3 <- ggplot(real_part3, aes(x=x, y=y)) +
geom_segment( aes(x=x, xend=x, y=0, yend=y) ) +
geom_point( size=3, color="blue", fill=alpha("lightblue", 0.1), alpha=0.7, shape=21, stroke=2) +
labs(title = "Estimation of x with lambda = 0.1", y = "Real part x_hat", x = "position") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(real_plot, real_plot1, real_plot2, real_plot3, ncol=1, nrow = 4)
As expected we fail to recover the true signal, because of the periodicity of the sampling method.
plot(error[,1], main = "Error for Lambda = 0.01", ylab = "Error", xlab = "Iterations")
plot(error[,2], main = "Error for Lambda = 0.05", ylab = "Error", xlab = "Iterations")
plot(error[,3], main = "Error for Lambda = 0.1", ylab = "Error", xlab = "Iterations")
The POCS algorithm fails to reconstruct the signal, since the sampling method of taking equidistant points periodises the signal.
im <- read.delim("brain_imag.txt", header=FALSE, sep=",")
im <- as.matrix(im)
real <- read.delim("brain_real.txt", header=FALSE, sep=",")
image(abs(im))
The temporal MRI signal (i.e. the raw data) directly samples the spatial frequency domain of the MR image. In other words, the measurements are the 2D-Fourier transform of the MR image i.e. a linear combination of the MR signal. These raw data are usually referred to as k-space. Compressive sensing MRI consists in sampling the k-space at a rate lower than the ShannonNyquist criterion. In other words, it consists in subsampling in the spatial-frequency domain. This is similar to what we studied in Part 2 and 3. The difference is that we studied 1D pure sparse signal and now we will study 2D not perfectly sparse signals.
First lets import the masks and pdfs :
mask_vardens <- read.delim("mask_vardens.txt", header=FALSE, sep=",")
mask_unif <- read.delim("mask_unif.txt", header=FALSE, sep=",")
pdf_vardens <- read.csv(file = "pdf_vardens.txt", header = F, sep = ',')
pdf_unif <- read.csv(file = "pdf_unif.txt", header = F, sep = ',')
M <- fft(im)/sqrt(length(im))
M_vardens <- as.matrix((M*mask_vardens)/pdf_vardens)
M_unif <- as.matrix((M*mask_unif)/pdf_unif)
im_vardens = fft(M_vardens, inverse=TRUE)/dim(M_vardens)[1]
im_unif = fft(M_unif, inverse=TRUE)/dim(M_unif)[1]
image(abs(im), main = "True image")
image(abs(im_vardens), main = "Vardens mask")
image(abs( (im_vardens - im)*10), main = "Error vardens mask")
image(abs(im_unif), main = "Unif mask")
image(abs((im_unif - im)*10), main = "Error unif mask")
Xhat <- M
i <- 1
while(i < 101) {
Xhat <- Xhat
xhat <- fft(Xhat, inverse = TRUE)/(dim(Xhat)[1])
for (i in 1:dim(xhat)[1]) {
xhat[i,] <- SoftThreshComplex(y = xhat[i,], lambda = 0.05)
}
Xhat <- (fft(xhat)/(dim(xhat)[1]))
Xhat <- Xhat*(M==0) + M
i <- i+1
}
xhat <- fft(Xhat, inverse = TRUE)/(dim(Xhat)[1])
image(abs(xhat))
SoftThreshComplex2D <- function(y, lambda) {
y <- as.vector(y)
xhat <- vector(length = length(y))
for( i in which((Mod(y)<lambda), arr.ind = TRUE)) {
xhat[i] <- 0
}
for( i in which((Mod(y)>lambda), arr.ind = TRUE) ) {
xhat[i] <- y[i]*(Mod(y[i])-lambda) / Mod(y[i])
}
matrix(xhat, ncol = 512, nrow = 512)
return(xhat)
}
Xhat <- M_vardens
i <- 1
while(i < 101) {
Xhat <- Xhat
xhat <- fft(Xhat, inverse = TRUE)/(dim(Xhat)[1])
xhat <- SoftThreshComplex2D(y = xhat, lambda = 0.01)
xhat <- matrix(xhat, ncol = 512, nrow = 512)
Xhat <- fft(xhat)/(dim(xhat)[1])
Xhat <- Xhat*(M_vardens==0) + M_vardens
i <- i+1
}
xhat <- fft(Xhat, inverse = TRUE)/dim(xhat)[1]
image(abs(xhat))
errorfun2d <- function(x, xhat){
return(sum(x^2-Im(xhat)^2))
}
Xhat <- M_vardens
i <- 1
error <- rep(0, 100)
while(i < 101) {
Xhat <- Xhat
xhat <- fft(Xhat, inverse = TRUE)/(dim(Xhat)[1])
for (k in 1:dim(xhat)[1]) {
xhat[k,] <- SoftThreshComplex(y = xhat[k,], lambda = 0.3)
}
Xhat <- (fft(xhat)/(dim(xhat)[1]))
Xhat <- Xhat*(M_vardens==0) + M_vardens
er <- rep(0, dim(xhat)[1])
for (j in 1:dim(xhat)[1]) {
er[j] <- errorfun2d(im[j,], xhat[j,])
}
error[i] <- sum(er)
i <- i+1
}
xhat <- fft(Xhat, inverse = TRUE)/(dim(Xhat)[1])
image(abs(xhat))
plot(error, main = "Error vardens", ylab = "Error", xlab = "Iterations")
Xhat <- M_unif
i <- 1
error <- rep(0, 100)
while(i < 101) {
Xhat <- Xhat
xhat <- fft(Xhat, inverse = TRUE)/(dim(Xhat)[1])
for (k in 1:dim(xhat)[1]) {
xhat[k,] <- SoftThreshComplex(y = xhat[k,], lambda = 0.25)
}
Xhat <- (fft(xhat)/(dim(xhat)[1]))
Xhat <- Xhat*(M_unif==0) + M_unif
er <- rep(0, dim(xhat)[1])
for (j in 1:dim(xhat)[1]) {
er[j] <- errorfun2d(im[j,], xhat[j,])
}
error[i] <- sum(er)
i <- i+1
}
xhat <- fft(Xhat, inverse = TRUE)/(dim(Xhat)[1])
image(abs(xhat))
plot(error, main = "Error unif", ylab = "Error", xlab = "Iterations")